/**
 * 预处理文件: data.txt  第1: int nSize(矩阵中非零个数)
 *N个int值表示每行中第一个数值的偏移 nSize个列 nSize个double类型的值
 *v1_b.txt 下面N个double类型的值
 *v1_x0.txt下面N个double类型的值
 *
**/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <omp.h>

#define Line 19
typedef struct
{
  int nx;
  int ny;
  int nz;
  int myrank;
}Info;
double *AV,*MV,*b,*x0,*tmp_x,*diag;
int *AJ;
int N;

Info InitInfo(int NX,int NY,int NZ,int PX,int PY,int PZ)
{
   Info result;
   int myrank;
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
   result.myrank=myrank;
   result.nx=NX;  result.ny=NY; result.nz=NZ;
   return result;
}

void LoadData(char *filename_A,char*filename_x0,char *filename_b)
{

    FILE *fp1=fopen(filename_A,"rb");
    if(fp1==NULL){
        printf("can't open read file\n");
        exit(0);
    }
    int size;
    fread(&size,4,1,fp1);
    AV=(double *)malloc(size* sizeof(double));
    AJ=(int*)malloc(size* sizeof(int));
    fread(AJ,4,size,fp1);
    fread(AV,8,size,fp1);
    fclose(fp1);

    FILE *fp2=fopen(filename_b,"rb");
    if(fp2==NULL){
        printf("can't open read file\n");
        exit(0);
    }
    b=(double *)malloc(N* sizeof(double));
    fread(b,8,N,fp2);
    fclose(fp2);

    FILE *fp3=fopen(filename_x0,"rb");
    if(fp3==NULL){
        printf("can't open read file\n");
        exit(0);
    }
    x0=(double *)malloc(N* sizeof(double));
    fread(x0,8,N,fp3);
    memcpy(tmp_x,x0,N*sizeof(double));
    fclose(fp3);
}

void getAmVector(double *result,double *data) //矩阵向量相乘
{
    int ind=0;
    int i,j;
    double tmp;
    #pragma omp parallel for private(i,j,tmp,ind) shared(AJ,data,AV,data,result)
    for(i=0;i<N;i++) //N-1
    {
        ind=i*Line;
        tmp=0.0;
        for(j=0;j<Line;j++){
            tmp=tmp+AV[ind+j]*data[AJ[ind+j]];
        }
        result[i]=tmp;
    }
}
double VectorDianJi(double *dataA,double *dataB)
{
    int i;
    double sum=0.0;
   #pragma omp parallel for private(i) shared(dataA,dataB) reduction(+:sum) 
   for(i=0;i<N;i++){
        sum=sum+dataA[i]*dataB[i];
    }
    return sum;
}

double checkSum(double *data)
{
    int ind=0;
    int i,j;
    double tmp;
    double sum=0.0;
    for(i=0;i<N;i++) //N-1行
    {
        ind=i*Line;
        tmp=0.0;
        for(j=0;j<Line;j++){
            tmp=tmp+AV[ind+j]*data[AJ[ind+j]];
        }
        sum=sum+(tmp-b[i])*(tmp-b[i]);
    }

    return sum;
}


void JacobiCal(double *Rm,double *R) //矩阵Rm=M^(-1)*R Jacobi迭代
{
    int ind=0;
    int i,j;
    double tmp,m_rst;
    for(i=0;i<N;i++){
        ind=i*Line;
        tmp=0.0;    m_rst=0.0;
        for(j=0;j<Line;j++){
            if(AJ[ind+j]==i){  //对角线上的元素
                tmp=AV[ind+j];
            }else{
                m_rst=m_rst+AV[ind+j]*R[AJ[ind+j]];
            }
        }
        Rm[i]=(R[i]-m_rst)/tmp;
    }
}

double getjk(int j,int k)//获取Ajk值
{
  
   int i,ind;
   ind=0;
   ind=j*Line;
   double reVal=0.0;
   for(i=0;i<Line;i++)
   {
      if(AJ[ind+i]==k){
         reVal=AV[ind+i];
         break;     
       }      
   }
   return reVal;
}
//ILU(0)
void getLUmatrixA()
{
    MV=(double *)malloc(N*Line*sizeof(double));
    int i,j,k;
    memcpy(MV,AV,N*Line*sizeof(double));    

    int ind;
    for(i=0;i<N;i++)
    {
       ind=i*Line;
       for(j=0;j<Line;j++){
          if(AJ[ind+j]==i){
             diag[i]=MV[ind+j];
             break;
          }
       }
       for(j=0;j<Line;j++)
       {
         if(AJ[ind+j]<i){
            MV[ind+j]=MV[ind+j]/diag[AJ[ind+j]];
            for(k=0;k<Line;k++){
               if(AJ[k+ind]>AJ[ind+j]){
                  MV[ind+k]=MV[ind+k]-MV[ind+j]*getjk(AJ[ind+j],AJ[ind+k]);
                  continue;
                }
            }
          }
       }
    }
}

void getMatrixLU(double *Rm,double *R)
{
   int i,j;
   double *z=(double*)malloc(N*sizeof(double));
   int ind;
   double tmp;
   for(i=0;i<N;i++)
   {
       ind=i*Line;
       for(j=0;j<Line;j++){
          if(AJ[ind+j]==i){
             diag[i]=MV[ind+j];
             break;
          }
       }
       tmp=0.0;
       for(j=0;j<Line;j++){
         if(AJ[ind+j]<i){
            tmp=tmp+MV[ind+j]*z[AJ[ind+j]];
         }
       }
      z[i]=(R[i]-tmp)/diag[i];
   }
   for(i=N-1;i>=0;i--)
   {
	ind=i*Line;
       tmp=0.0;
       for(j=0;j<Line;j++){
         if(AJ[ind+j]>i){
            tmp=tmp+MV[ind+j]*Rm[AJ[ind+j]];
         }
       }
       Rm[i]=z[i]-tmp/diag[i];
   }
   free(z);
}

//D-ILU(0)分解
void factorDLU()
{
   int ind,i,j;
   ind=0;
   for(i=0;i<N;i++){
       ind=i*Line;
     for(j=0;j<Line;j++){
	if(AJ[ind+j]==i)
	    diag[i]=AV[ind+j];
	    break;
	}
    }    

    ind=0;
    for(i=0;i<N;i++)
    {
       diag[i]=1.0/diag[i];
       ind=i*Line;
       for(j=0;j<Line;j++){
         if(AJ[ind+j]>i){
           double tmp=getjk(AJ[ind+j],i);
	   if(tmp!=0){
	      diag[AJ[ind+j]]=diag[AJ[ind+j]]-tmp*diag[i]*AV[ind+j];	  	
           }    
	}      
      }
   }
}
void getDLUCal(double *Rm,double *R)
{
   int i,j,ind;
   double *z=(double*)malloc(N*sizeof(double)); 
   double tmp;
   ind=0;
   for(i=0;i<N;i++)
   {
       ind=i*Line;
       tmp=0.0;
       for(j=0;j<Line;j++){
         if(AJ[ind+j]<i){
            tmp=tmp+AV[ind+j]*z[AJ[ind+j]];
         }
       }
      z[i]=R[i]-tmp*diag[i];
   }
   ind=0;
   for(i=N-1;i>=0;i--)
   {
       ind=i*Line;
       tmp=0.0;
       for(j=0;j<Line;j++){
         if(AJ[ind+j]>i){
            tmp=tmp+AV[ind+j]*Rm[AJ[ind+j]];
         }
       }
       Rm[i]=z[i]-tmp*diag[i];
   }
   free(z); 
}

void calGCR(int flag)
{
     memcpy(x0,tmp_x,N*sizeof(double));
    int k=5;
    int i,j;
    double *R=(double*)malloc(N*sizeof(double));
    double *Rm=(double*)malloc(N*sizeof(double));
    double *p=(double*)malloc(N*sizeof(double));
    double *Ap=(double*)malloc(N* sizeof(double));
    double *ARm=(double*)malloc(N* sizeof(double));
    double *Apj=(double*)malloc(N* sizeof(double));

    double *pj[5];
    double betaij[5];
    for(i=0;i<5;i++){
        pj[i]=(double*)malloc(N*sizeof(double));
        betaij[i]=0;
    }
    getAmVector(R,x0);
    for(i=0;i<N;i++){   //get the R=b-Ax0
        R[i]=b[i]-R[i];
    }
   if(flag==0) memcpy(Rm,R,N*sizeof(double));
   if(flag==1) JacobiCal(Rm,R);
   if(flag==2) getDLUCal(Rm,R);
   if(flag==3) getMatrixLU(Rm,R);    

    memcpy(p,Rm,N*sizeof(double));
    double ai,checkVal;
    int ind=0,m=0;
    i=1;
    do{
        getAmVector(Ap,p);   //get Api-1 向量 
        
        for(j=0;j<N;j++){
          pj[(i-1)%k][j]=p[j];      //存储5个子空间向量 
        }

        ai=VectorDianJi(R,Ap)/VectorDianJi(Ap,Ap);
        for(j=0;j<N;j++){       //xi=x(i-1)+a(i-1)*p(i-1)
            x0[j]=x0[j]+ai*p[j];
            R[j]=R[j]-ai*Ap[j];
       }
	 if(flag==0) memcpy(Rm,R,N*sizeof(double));
  	 if(flag==1) JacobiCal(Rm,R);
   	 if(flag==2) getDLUCal(Rm,R);
         if(flag==3) getMatrixLU(Rm,R);
     //   JacobiCal(Rm,R);        //预处理
    //    getDLUCal(Rm,R);
   //     getMatrixLU(Rm,R);    
        ind=0;
        for(j=((i-1)/k)*k;j<i;j++){
            getAmVector(ARm,Rm);
            getAmVector(Apj,pj[ind]);
            betaij[ind]=-VectorDianJi(ARm,Apj)/VectorDianJi(Apj,Apj);
            ind++;
        }

       for(j=0;j<N;j++){
          p[j]=Rm[j];
       }     
        ind=0;
        for(j=((i-1)/k)*k;j<i;j++){
            for(m=0;m<N;m++){
                p[m]=p[m]+betaij[ind]*pj[ind][m];
            }    
            ind++;
        }
       checkVal=checkSum(x0);
     //  printf("checkSum---->%.30lf\n",checkVal); 
       i++;
    }while(checkVal>0.0000000001);
    printf("The final checkSum---->%.30lf\n",checkVal); 
    printf("Iterator times------>%d\n",(i-1));
    free(R);free(Rm);free(p);free(Ap);free(ARm); free(Apj);
    for(i=0;i<5;i++){
       free(pj[i]);
    }
}

int main(int argc, char **argv) {
    AJ=NULL; AV=NULL; b=NULL; x0=NULL;

    MPI_Init(&argc,&argv);
    int myrank;
    int NX,NY,NZ;
    int PX,PY,PZ;
    NX=atoi(argv[1]);NY = atoi(argv[2]); NZ = atoi(argv[3]); 
    PX=atoi(argv[4]);PY = atoi(argv[5]); PZ = atoi(argv[6]); 
    char *filename_A=argv[7];
    char *filename_x0=argv[8];
    char *filename_b=argv[9];
    if(PZ!=1) return 1;
    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
    Info info=InitInfo(NX,NY,NZ,PX,PY,PZ); 
   if(myrank==0){
  
     N=NX*NY*NZ;
     tmp_x=(double*)malloc(N*sizeof(double));
     diag=(double*)malloc(N*sizeof(double));
     LoadData(filename_A,filename_x0,filename_b);
     printf("the Init checkSum---->%.30lf\n", checkSum(x0));
     // 0: 无预处理 1表示Jacobi预处理 2表示D-ILU(0)预处理 3表示ILU(0)预处理
     calGCR(0);
     calGCR(1);     
     factorDLU();
     calGCR(2);
     getLUmatrixA();
     calGCR(3);
    free(tmp_x);
    free(AJ);free(AV);free(b);free(x0);
    free(diag);
  }
   MPI_Finalize();
   return 0;
}

































